Descenso de Gradiente y Auto-diferenciación

Ahora queremos estudiar cómo encontrar los parámetros óptimos de un modelo para minimizar una función de pérdida. Existen muchos algoritmos de optimización. Uno de los que converge en menos pasos es el método de Newton, pero su costo computacional crece rápidamente con el número de parámetros. Las redes neuronales pueden involucrar billones de parámetros, lo que hace que el método de Newton y otros similares sean imposibles de aplicar. Por lo tanto se usa el Descenso de Gradiente (Gradient Descent, GD) al entrenar redes neuronales.

El descenso de gradiente consiste en tomar el gradiente de la función de pérdida en el espacio de parámetros. La idea es ir en la dirección opuesta al gradiente, ya que el gradiente apunta en la dirección de máximo crecimiento de la función (y su opuesto en la dirección de máximo descenso). Cuando llegamos a un mínimo de la función, el gradiente será cero.

En esta sección veremos primero cómo aplicar el descenso de gradiente para un perceptrón multicapa asumiendo que el gradiente de la función de pérdida con respecto a los parámetros de la red es conocido. Luego veremos cómo calcular el gradiente de la función de pérdida con respecto a los parámetros de la red de manera eficiente. Finalmente discutiremos algunos detalles de cómo escoger el punto inicial en el espacio de parámetros para inciar el descenso de gradiente.

Descenso de Gradiente

Ejemplo de Regresión Lineal Simple con Descenso de Gradiente

Antes de aplicar el descenso de gradiente a una red neuronal, empezamos por hacerlo en un caso sencillo que es fácil de visualizar: la regresión lineal simple. En este caso queremos ajustar la pendiente y el intercepto de una recta \(y = m x + b\) a un conjunto de datos \((x_i, y_i)\). La función de pérdida (MSE) a minimizar es: \[ \mathcal{J}(m,b) = \frac{1}{N}\sum_{i=1}^N (y_i - (m x_i + b))^2. \] donde \(N\) es el número total de puntos de datos. Notemos que la función de pérdida es una función de dos variables, \(m\) y \(b\). Por eso la podemos graficar fácilmente.

Código 
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Generar muestra aleatoria
np.random.seed(42)
N = 100
x = np.random.uniform(-10, 10, N)
m_verdadero, b_verdadero = 2, -1
y = m_verdadero * x + b_verdadero + np.random.normal(0, 2, N)

# Definir función de pérdida
def perdida_ecm(m, b, x, y):
    return np.mean((y - (m * x + b))**2)

# Crear malla para gráfico 3D y heatmap
rango_m = np.linspace(-3, 7, 100)
rango_b = np.linspace(-24, 24, 100)
M, B = np.meshgrid(rango_m, rango_b)
Z = np.array([perdida_ecm(m, b, x, y) for m, b in zip(M.ravel(), B.ravel())]).reshape(M.shape)

# Configurar el layout de los subplots
fig = plt.figure(figsize=(6, 12))

# Graficar superficie 3D
ax1 = fig.add_subplot(211, projection='3d')
superficie = ax1.plot_surface(M, B, Z, cmap='viridis')
ax1.set_xlabel('m (pendiente)')
ax1.set_ylabel('b (intercepto)')
ax1.set_zlabel('Pérdida ECM')
ax1.set_title('Función de Pérdida 3D')

# Graficar heatmap 2D
ax2 = fig.add_subplot(212)
heatmap = ax2.imshow(Z, extent=[rango_m.min(), rango_m.max(), rango_b.min(), rango_b.max()], 
                     origin='lower', cmap='viridis', aspect='auto')
contours = ax2.contour(M, B, Z, colors='black', levels=10)
ax2.set_xlabel('m (pendiente)')
ax2.set_ylabel('b (intercepto)')
ax2.set_title('Heatmap de la Función de Pérdida')

plt.tight_layout()
plt.show()
Figura 1: Función de pérdida para regresión lineal

En la Figura 1 se muestra la función de pérdida para regresión lineal. Vemos que esta es mínima en el punto que corresponde a los valores verdaderos de \(m\) y \(b\) y que es convexa, es decir, que cualquier línea recta que conecte dos puntos está por encima de la superficie.

Utilizamos el Descenso de Gradiente (Gradient Descent, GD) para encontrar los valores óptimos de \(m\) y \(b\). La actualización en cada iteración se realiza siguiendo la dirección opuesta al gradiente de la función de pérdida:

\[ \begin{aligned} m &\leftarrow m - \eta \,\frac{\partial \mathcal{L}}{\partial m},\\ b &\leftarrow b - \eta \,\frac{\partial \mathcal{L}}{\partial b}, \end{aligned} \]

Aquí, \(\eta\) es la tasa de aprendizaje (learning rate), un hiperparámetro que controla el tamaño de los pasos que damos en cada actualización. Las derivadas parciales (el gradiente) indican la dirección de máximo ascenso de la pérdida, y su inverso es la dirección de máximo descenso. Al restar el gradiente multiplicado por \(\eta\), nos movemos un poco hacia el mínimo. En este caso podemos calcular fácilmente el gradiente de la función de pérdida con respecto a \(m\) y \(b\): \[ \frac{\partial \mathcal{J}}{\partial m} = -\frac{2}{N}\sum_{i=1}^N (y_i - (m x_i + b))\,x_i, \] \[ \frac{\partial \mathcal{J}}{\partial b} = -\frac{2}{N}\sum_{i=1}^N (y_i - (m x_i + b)). \]

Código 
import matplotlib.animation as animation
from IPython.display import HTML

# Definir función de pérdida y sus derivadas
def perdida_ecm(m, b, x, y):
    return np.mean((y - (m * x + b))**2)

def derivada_m(m, b, x, y):
    return -2 * np.mean((y - (m * x + b)) * x)

def derivada_b(m, b, x, y):
    return -2 * np.mean(y - (m * x + b))

# Crear malla para gráfico 3D y heatmap
rango_m = np.linspace(0, 4, 100)
rango_b = np.linspace(-5.5, 2, 100)
M, B = np.meshgrid(rango_m, rango_b)
Z = np.array([perdida_ecm(m, b, x, y) for m, b in zip(M.ravel(), B.ravel())]).reshape(M.shape)

# Configurar figura
fig, ax = plt.subplots(figsize=(10, 6))
heatmap = ax.imshow(Z, extent=[rango_m.min(), rango_m.max(), rango_b.min(), rango_b.max()], 
                    origin='lower', cmap='viridis', aspect='auto')
contours = ax.contour(M, B, Z, colors='black', levels=10)
ax.set_xlabel('m (pendiente)')
ax.set_ylabel('b (intercepto)')
ax.set_title('Descenso de Gradiente')

# Ejecutar el descenso de gradiente
m_inicio, b_inicio = 3, -2.5  # Punto inicial
eta = 0.015  # Tasa de aprendizaje
pasos = 15  # Número de pasos a mostrar

# Calcular la trayectoria completa de antemano
trayectoria_m = [m_inicio]
trayectoria_b = [b_inicio]

m_actual, b_actual = m_inicio, b_inicio
for _ in range(pasos):
    grad_m = derivada_m(m_actual, b_actual, x, y)
    grad_b = derivada_b(m_actual, b_actual, x, y)
    
    # Actualizar parámetros
    m_actual -= eta * grad_m
    b_actual -= eta * grad_b
    
    # Guardar trayectoria
    trayectoria_m.append(m_actual)
    trayectoria_b.append(b_actual)

# Inicializar elementos de la animación
line, = ax.plot([], [], 'ro-', linewidth=2, markersize=8)
points = ax.scatter([], [], color='red', s=50)

# Función de inicialización
def init():
    line.set_data([], [])
    points.set_offsets(np.empty((0, 2)))
    return line, points

# Función de actualización para cada frame
def update(frame):
    line.set_data(trayectoria_m[:frame+1], trayectoria_b[:frame+1])
    points.set_offsets(np.column_stack((trayectoria_m[:frame+1], trayectoria_b[:frame+1])))
    return line, points

# Crear animación
ani = animation.FuncAnimation(fig, update, frames=len(trayectoria_m),
                              init_func=init, blit=True, interval=300)

HTML(ani.to_jshtml())
(a) Descenso de gradiente para regresión lineal
(b)
Figura 2

En la Figura 2 se muestra el descenso de gradiente para regresión lineal. Graficamos unos pocos pasos del descenso de gradiente. Ya podemos apreciar un posible problema: El descenso avanza rápidamente en la dirección de la pendiente, pero luego se detiene y avanza lentamente en la dirección del intercepto. Esto es porque el gradiente es más pronunciado en la dirección de la pendiente que en la dirección del intercepto. Veremos más adelante cómo intentar mejorar esto.

Ejemplo 2D de modelo Gabor: paisajes con mínimos locales

El ejemplo de arriba fue para una función de pérdida convexa. En general las funciones de pérdida de las redes neuronales no tienen esta propiedad. Por lo tanto, el “paisaje” de la pérdida puede ser complejo, presentando múltiples mínimos locales y puntos de silla (saddle points). Esto dificulta la optimización, ya que el Descenso de Gradiente puede quedar atrapado en un mínimo local que no es el mínimo global. Usaremos una función sencilla con esas características llamada función de Gabor para ilustrar el problema.

Usamos una función Gabor 1D como modelo, donde ajustaremos dos de sus parámetros. Consideramos \(\sigma\) y \(\omega\) como parámetros fijos conocidos. Optimizaremos la amplitud \(A\) y la fase \(\phi\): \[ f(x; \phi_0,\phi_1) = e^{-\frac{(\phi_0 + \phi_1 x)^2}{32}}\cos(\phi_0 + \phi_1 x + \phi). \] Ajustamos \((\phi_0,\phi_1)\) a datos \((x_i, y_i)\) generados usando valores verdaderos \((\phi_0^*,\phi_1^*)\), en este ejemplo.

Minimizamos el Error Cuadrático Medio entre las predicciones del modelo Gabor y los datos \(y_i\). Siendo \(N\) el número de puntos \(x_i\): \[ \mathcal{L}(\phi_0,\phi_1) = \frac{1}{N}\sum_{i=1}^N (y_i - f(x_i; \phi_0,\phi_1))^2. \] La naturaleza oscilatoria del término \(\cos(\phi_0 + \phi_1 x + \phi)\), multiplicada por la amplitud \(A\), crea un paisaje de pérdida no convexo con múltiples valles (mínimos locales) y puntos de silla.

Para entender el paisaje de la pérdida, generamos un gráfico de contorno de \(\mathcal{L}(\phi_0,\phi_1)\) evaluando la pérdida en una rejilla de valores de \(\phi_0\) y \(\phi_1\). Este gráfico muestra las regiones de baja pérdida (valles) y las zonas de puntos de silla. Adicionalmente, se pueden trazar las trayectorias seguidas por el algoritmo de Descenso de Gradiente partiendo desde diferentes puntos iniciales \((\phi_{0,0}, \phi_{1,0})\).

Código 
# Definir función de Gabor y generar datos
def gabor(x, phi_0, phi_1, sigma=5.0):
    return np.exp(-(phi_0 + phi_1*x)**2/(2*sigma**2)) * np.sin(phi_0 + phi_1*x)

def d_gabor_dx(x, sigma=5.0):
    return - x * np.exp(-x**2/(2*sigma**2))/sigma**2 * np.sin(x) + \
        np.exp(-x**2/(2*sigma**2)) * np.cos(x)

# Generar datos sintéticos
np.random.seed(42)
N = 100
x = np.linspace(-3, 3, N)
phi_0_verdadero, phi_1_verdadero = 2.0, 1.0
sigma = 5.0  # Parámetros fijos
y = gabor(x, phi_0_verdadero, phi_1_verdadero, sigma) + np.random.normal(0, 1, N)

# Definir función de pérdida y sus derivadas
def perdida_ecm(phi_0, phi_1, x, y, sigma=5.0):
    y_pred = gabor(x, phi_0, phi_1, sigma)
    return np.mean((y - y_pred)**2)

def derivada_phi_0(phi_0, phi_1, x, y, sigma=5.0):
    y_pred = gabor(x, phi_0, phi_1, sigma)
    return -2*np.mean((y - y_pred) * d_gabor_dx(phi_0 + phi_1*x, sigma))

def derivada_phi_1(phi_0, phi_1, x, y, sigma=5.0):
    y_pred = gabor(x, phi_0, phi_1, sigma)
    return -2*np.mean((y - y_pred) * d_gabor_dx(phi_0 + phi_1*x, sigma)*x)

# Crear malla para gráfico de contorno
rango_phi_0 = np.linspace(-1, 4, 100)
rango_phi_1 = np.linspace(-2, 4, 100)
Phi_0_grid, Phi_1_grid = np.meshgrid(rango_phi_0, rango_phi_1)
Z = np.array([perdida_ecm(p_0, p_1, x, y) for p_0, p_1 in zip(Phi_0_grid.ravel(), Phi_1_grid.ravel())]).reshape(Phi_0_grid.shape)

# Configurar figura
fig, ax = plt.subplots(figsize=(10, 6))
heatmap = ax.imshow(Z, extent=[rango_phi_0.min(), rango_phi_0.max(), rango_phi_1.min(), rango_phi_1.max()], 
                    origin='lower', cmap='viridis', aspect='auto')
contours = ax.contour(Phi_0_grid, Phi_1_grid, Z, colors='black', levels=10)
ax.set_xlabel(r'$\phi_0$')
ax.set_ylabel(r'$\phi_1$')
ax.set_title('Descenso de Gradiente en Modelo Gabor')

# Definir puntos iniciales para diferentes trayectorias
puntos_iniciales = [
    (2.3, 1.0),    # (φ_0, φ_1) - Punto 1
    (1.0, 3.0),    # Punto 2
    (3.0, 0.0),    # Punto 3
    (1.0, -1.0),   # Punto 4
    (2.5, -1.5)    # Punto 5
]

colores = ['r', 'b', 'g', 'y', 'c']
eta = 0.1  # Tasa de aprendizaje
pasos = 30  # Número de pasos a mostrar

# Calcular todas las trayectorias de antemano
todas_trayectorias = []

for phi_0_inicio, phi_1_inicio in puntos_iniciales:
    trayectoria_phi_0 = [phi_0_inicio]
    trayectoria_phi_1 = [phi_1_inicio]
    
    phi_0_actual, phi_1_actual = phi_0_inicio, phi_1_inicio
    for _ in range(pasos):
        grad_phi_0 = derivada_phi_0(phi_0_actual, phi_1_actual, x, y)
        grad_phi_1 = derivada_phi_1(phi_0_actual, phi_1_actual, x, y)
        
        # Actualizar parámetros
        phi_0_actual -= eta * grad_phi_0
        phi_1_actual -= eta * grad_phi_1
        
        # Guardar trayectoria
        trayectoria_phi_0.append(phi_0_actual)
        trayectoria_phi_1.append(phi_1_actual)
    
    todas_trayectorias.append((trayectoria_phi_0, trayectoria_phi_1))

# Inicializar elementos de la animación
lines = []
points_collection = []

for i, color in enumerate(colores):
    line, = ax.plot([], [], f'{color}-', linewidth=2, markersize=8)
    points = ax.scatter([], [], color=color, s=50)
    lines.append(line)
    points_collection.append(points)

# Función de inicialización
def init():
    for line in lines:
        line.set_data([], [])
    for points in points_collection:
        points.set_offsets(np.empty((0, 2)))
    return lines + points_collection

# Función de actualización para cada frame
def update(frame):
    for i, ((trayectoria_phi_0, trayectoria_phi_1), line, points) in enumerate(zip(todas_trayectorias, lines, points_collection)):
        line.set_data(trayectoria_phi_0[:frame+1], trayectoria_phi_1[:frame+1])
        points.set_offsets(np.column_stack((trayectoria_phi_0[:frame+1], trayectoria_phi_1[:frame+1])))
    
    return lines + points_collection

# Crear animación
ani = animation.FuncAnimation(fig, update, frames=pasos+1,
                              init_func=init, blit=True, interval=300)

HTML(ani.to_jshtml())
(a) Descenso de gradiente para modelo Gabor
(b)
Figura 3

En la Figura 3 se muestra el paisaje de la pérdida para el modelo Gabor, y las trayectorias de convergencia para diferentes puntos iniciales. Vemos que diferentes puntos van hacia diferentes mínimos locales, lo que dificulta la convergencia al mínimo global.

Descenso de Gradiente Estocástico (SGD) y por mini-lotes

Cuando el conjunto de datos es muy grande, calcular el gradiente sobre todos los \(N\) ejemplos en cada paso puede ser computacionalmente muy costoso. Por ello, se utilizan variantes del GD que usan un subconjunto de los datos para aproximar el gradiente a cada paso. Esta será sólo una aproximación, y nos da una estima ruidosa del gradiente de la pérdida ideal. Pero esto tiene una gran ventaja: Ese ruido puede empujar el algoritmo a saltar de mínimos locales poco profundos. Veamos dos variaciones de esa idea:

  • Descenso de Gradiente Estocástico (SGD) Realiza una actualización de parámetros por cada ejemplo individual \((x_i, y_i)\). En cada paso, el gradiente se estima usando solo ese único ejemplo: \(\theta \leftarrow \theta - \eta \nabla_{\theta} \mathcal{L}_{\text{indiv}}(\theta; x_i, y_i)\), donde \(\mathcal{L}_{\text{indiv}}\) es la pérdida calculada para ese punto. Luego repite para todos los puntos \((x_i, y_i)\) y vuelve a empezar.
    • Ventajas: Actualizaciones muy rápidas y bajo coste computacional por actualización. La alta varianza en la estimación del gradiente puede ayudar a escapar de mínimos locales “malos”.
    • Desventajas: Al usar un único dato para calcular el gradiente, la trayectoria de convergencia es muy ruidosa (alta varianza). Puede que nunca converja exactamente al mínimo, sino que oscile alrededor de este.
  • Descenso de Gradiente por Mini-lotes: Es un compromiso entre los dos anteriores y el método más común en Deep Learning. El gradiente se calcula sobre un pequeño subconjunto aleatorio de datos llamado mini-batch (lote), de tamaño \(B\) (e.g., \(B=32, 64, 128\)). Luego repite para todos los lotes y vuelve a empezar.
    • Ventajas: Reduce significativamente la varianza del gradiente comparado con SGD, llevando a una convergencia más estable. Mantiene una alta eficiencia computacional aprovechando la paralelización de hardware (GPUs).
    • Desventajas: Introduce un nuevo hiperparámetro (tamaño del batch, \(B\)).

En la práctica, el descenso de gradiente por mini-lotes es el método más común y eficiente. Es el usado para entrenar las redes neuronales famosas que han impulsado el interés en el aprendizaje automático.

A continuación el pseudocódigo para el descenso de gradiente por mini-lotes:

# Inicializar parámetros
theta = theta_0

# Bucle de entrenamiento
for epoch in range(num_epochs):
    # Mezclar datos
    np.random.shuffle(data)
    
    # Bucle sobre lotes
    for batch in data:
        # Calcular gradiente
        grad = compute_gradient(batch)
        
        # Actualizar parámetros
        theta = theta - eta * grad

Momentum

Otra estrategia para atacar el problema de los mínimos locales y la oscilación alrededor de mínimos, es cambiar el proceso de aprendizaje. Es decir, modificamos la forma como se actualizan los pesos a cada paso.

El método de Momentum introduce un término de “velocidad” \(v\) que acumula una media móvil exponencial de los gradientes pasados. Ayuda a acelerar el descenso en direcciones consistentes y amortigua las oscilaciones en direcciones donde el gradiente cambia rápidamente (típico en valles estrechos y alargados).

La actualización se modifica como sigue (donde \(\theta\) representa los parámetros del modelo, como \(m\) y \(b\)):

\[ \begin{aligned} v_t &\leftarrow \beta\,v_{t-1} + (1-\beta)\,\nabla\mathcal{L}(\theta_{t-1}),\\ \theta_t &\leftarrow \theta_{t-1} - \eta\,v_t. \end{aligned} \]

  • \(v_t\) es el vector de velocidad en el paso \(t\).
  • \(\beta\) es el coeficiente de momentum (usualmente cercano a 0.9). Controla cuánto del momentum anterior se conserva. Un \(\beta\) alto significa que los gradientes pasados tienen más influencia.
  • \(\nabla\mathcal{L}(\theta_{t-1})\) es el gradiente calculado en el paso actual (puede ser batch, mini-batch o stochastic).

A cantidades que se acumulan de la manera como lo hace \(v_t\) se les llama media móvil exponencial.

Esto es burdamente similar a una bola pesada rodando por una pendiente. Acumula momento, lo que le permite pasar por pequeños baches y acelerar en las bajadas consistentes. A pesar de llamarse momentum, no es exactamente el momentum definido en la física.

A continuación el pseudocódigo para el descenso de gradiente con momentum en el problema de regresión lineal:

# Inicializar velocidades
v_m, v_b = 0.0, 0.0
beta = 0.9 # Coeficiente de momentum
eta = 0.01 # Tasa de aprendizaje

# Bucle de entrenamiento (puede ser sobre épocas y mini-batches)
for iteration in range(num_iterations):
    # Calcular gradientes (usando batch, mini-batch, o SGD)
    # grad_m = ∂L/∂m, grad_b = ∂L/∂b
    grad_m, grad_b = compute_gradient(m, b, data_batch)

    # Actualizar velocidades
    v_m = beta * v_m + (1 - beta) * grad_m
    v_b = beta * v_b + (1 - beta) * grad_b

    # Actualizar parámetros usando las velocidades
    m -= eta * v_m
    b -= eta * v_b

Adam (Adaptive Moment Estimation)

Cuando hicimos el ejemplo de regresión lineal, vimos que la convergencia hacia el mínimo en una dimensión era mucho más rápida que en la otra. Adam es un optimizador que intenta resolver este problema ajustando automáticamente la tasa de aprendizaje para cada parámetro individualmente.

Adam es un optimizador muy popular que combina dos ideas principales:

  1. Momentum: Utiliza una media móvil exponencial del gradiente (primer momento, \(m_t\)), similar al método de Momentum.
  2. Adaptación de Tasa de Aprendizaje: Utiliza una media móvil exponencial del cuadrado del gradiente (segundo momento, \(v_t\)). Esto reduce la diferencia entre las tasas de aprendizaje más rápidas y las más lentas.

Además, aplica una corrección del sesgo (bias correction) a ambas estimaciones de momentos para contrarrestar su tendencia a estar sesgadas hacia cero durante las primeras iteraciones.

Las ecuaciones de actualización son: \[ \begin{aligned} m_t &= \beta_1 m_{t-1} + (1-\beta_1)\,g_t & \text{(Estimación del 1er momento)} \\ v_t &= \beta_2 v_{t-1} + (1-\beta_2)\,g_t^2 & \text{(Estimación del 2do momento)} \\ \hat m_t &= \frac{m_t}{1-\beta_1^t} & \text{(Corrección del sesgo 1er momento)} \\ \hat v_t &= \frac{v_t}{1-\beta_2^t} & \text{(Corrección del sesgo 2do momento)} \\ \theta_t &\leftarrow \theta_{t-1} - \eta\,\frac{\hat m_t}{\sqrt{\hat v_t}+\epsilon} & \text{(Actualización del parámetro } \theta \text{)} \end{aligned} \]

Donde: * \(t\) es el número de iteración (paso de actualización). * \(g_t = \nabla_{\theta} \mathcal{L}(\theta_{t-1})\) es el gradiente en el paso \(t\). * \(\beta_1\) es el coeficiente de decaimiento para el primer momento (usualmente ~0.9). * \(\beta_2\) es el coeficiente de decaimiento para el segundo momento (usualmente ~0.999). * \(\eta\) es la tasa de aprendizaje global (e.g., 1e-3). * \(\epsilon\) es una constante pequeña (usualmente ~1e-8) para evitar la división por cero y mejorar la estabilidad numérica.

A continuación se muestra un pseudocódigo para el descenso de gradiente con Adam en el problema de regresión lineal.

# Inicializar momentos
m_m, m_b = 0.0, 0.0  # Primer momento (media móvil del gradiente)
v_m, v_b = 0.0, 0.0  # Segundo momento (media móvil del cuadrado del gradiente)
beta1 = 0.9          # Coeficiente de decaimiento para el primer momento
beta2 = 0.999        # Coeficiente de decaimiento para el segundo momento
epsilon = 1e-8       # Constante para evitar división por cero
eta = 0.01           # Tasa de aprendizaje
t = 0                # Contador de iteraciones

# Bucle de entrenamiento
for iteration in range(num_iterations):
    t += 1
    
    # Calcular gradientes (usando batch, mini-batch, o SGD)
    grad_m, grad_b = compute_gradient(m, b, data_batch)
    
    # Actualizar primer momento (media móvil del gradiente)
    m_m = beta1 * m_m + (1 - beta1) * grad_m
    m_b = beta1 * m_b + (1 - beta1) * grad_b
    
    # Actualizar segundo momento (media móvil del cuadrado del gradiente)
    v_m = beta2 * v_m + (1 - beta2) * (grad_m**2)
    v_b = beta2 * v_b + (1 - beta2) * (grad_b**2)
    
    # Corregir el sesgo
    m_m_corrected = m_m / (1 - beta1**t)
    m_b_corrected = m_b / (1 - beta1**t)
    v_m_corrected = v_m / (1 - beta2**t)
    v_b_corrected = v_b / (1 - beta2**t)
    
    # Actualizar parámetros
    m -= eta * m_m_corrected / (np.sqrt(v_m_corrected) + epsilon)
    b -= eta * m_b_corrected / (np.sqrt(v_b_corrected) + epsilon)

Diferenciación Automática (AD)

Vimos que el descenso de gradiente y los algoritmos inspirados en este usan la derivada de la función de pérdida con respecto a los parámetros de la red. Las redes neuronales tienen hasta billones de parámetros, por lo que calcular la derivada de la función de pérdida con respecto a cada uno de ellos puede ser bastante costoso. Sin embargo, se han desarrollado técnicas que atacan este problema de forma bastante ingeniosa.

La diferenciación automática (AD) es un conjunto de técnicas que permite calcular de forma exacta y eficiente derivadas de funciones definidas por programas informáticos. Es crucial distinguirla de la diferenciación simbólica (que puede llevar a expresiones muy complejas) y de la diferenciación numérica (que usa aproximaciones por diferencias finitas y sufre de errores de truncamiento y redondeo). AD calcula derivadas exactas (hasta la precisión de la máquina) con una eficiencia comparable a la evaluación de la propia función.

Algoritmo de Backpropagation (Modo Reverso de AD)

En redes neuronales y otros modelos con muchos parámetros, la diferenciación automática se implementa de forma eficiente mediante el modo reverso (reverse mode), cuyo algoritmo más conocido es backpropagation. Existe también el modo directo (forward mode), que es menos eficiente pero puede ser útil en ciertas situaciones, y no lo veremos aquí por falta de tiempo.

El modo reverso calcula todas las derivadas de una única salida (escalar, como la función de pérdida \(\mathcal{J}\)) con respecto a todas las entradas y parámetros intermedios. Lo hace propagando las derivadas hacia atrás, desde la salida hacia la entrada, aplicando eficientemente la regla de la cadena sobre el grafo computacional definido por la evaluación de la función (el forward pass).

El proceso tiene dos fases:

  1. Forward pass (Paso hacia adelante)
    • Partimos de la entrada \(x\) y calculamos, capa a capa, las activaciones intermedias y la salida final (y la pérdida \(\mathcal{J}\)). \[ z^{(1)} = W^{(1)} x + b^{(1)}, \quad a^{(1)} = \sigma(z^{(1)}), \quad z^{(2)} = W^{(2)} a^{(1)} + b^{(2)}, \dots, \quad \mathcal{J} = \text{Loss}(\text{output}, \text{target}) \]
    • Se almacenan en memoria todas las variables intermedias (\(z^{(\ell)}\), \(a^{(\ell)}\), \(W^{(\ell)}\), \(x\)) que influyen en el resultado final, ya que serán necesarias para calcular las derivadas en el paso inverso.
  2. Backward pass (Paso hacia atrás)
    • Se inicializa con el gradiente de la pérdida respecto a sí misma (\(\frac{\partial \mathcal{J}}{\partial \mathcal{J}} = 1\)) o respecto a la salida final de la red antes de la función de pérdida.
    • Se propaga el gradiente hacia atrás capa por capa (\(\ell = L, L-1, \dots, 1\)). Para cada capa \(\ell\):
      • Se parte del gradiente respecto a la activación de la capa, \(\frac{\partial \mathcal{J}}{\partial a^{(\ell)}}\) (obtenido de la capa \(\ell+1\), o inicial si \(\ell=L\)).
      • Se calcula el gradiente respecto a la entrada lineal de la capa \(\ell\), denotado \(\delta^{(\ell)}\): \[ \delta^{(\ell)} = \frac{\partial \mathcal{J}}{\partial z^{(\ell)}} = \frac{\partial \mathcal{J}}{\partial a^{(\ell)}} \odot \sigma'(z^{(\ell)}) \] donde \(\odot\) es el producto elemento a elemento (Hadamard) y \(\sigma'\) es la derivada de la función de activación \(\sigma\) aplicada elemento a elemento a \(z^{(\ell)}\). En el caso de ReLU, \(\sigma'(z^{(\ell)}) = 1\) para \(z^{(\ell)} > 0\) y \(\sigma'(z^{(\ell)}) = 0\) para \(z^{(\ell)} \leq 0\).
      • Se calculan los gradientes respecto a los parámetros de la capa \(\ell\) usando \(\delta^{(\ell)}\) y la activación de la capa anterior \(a^{(\ell-1)}\) (guardada durante el forward pass): \[ \frac{\partial \mathcal{J}}{\partial W^{(\ell)}} = \delta^{(\ell)} \, (a^{(\ell-1)})^T \] \[ \frac{\partial \mathcal{J}}{\partial b^{(\ell)}} = \delta^{(\ell)} \quad \text{(sumado sobre las muestras del batch si aplica)} \]
      • Se calcula el gradiente que se propagará a la capa anterior \(\ell-1\): \[ \frac{\partial \mathcal{J}}{\partial a^{(\ell-1)}} = (W^{(\ell)})^T \, \delta^{(\ell)} \] Este \(\frac{\partial \mathcal{J}}{\partial a^{(\ell-1)}}\) será el punto de partida para los cálculos de la capa \(\ell-1\).
    • Este proceso continúa hasta calcular los gradientes respecto a todos los parámetros \(W^{(\ell)}, b^{(\ell)}\) de la red (y respecto a la entrada \(x\) si fuera necesario).

Eficiencia: Backpropagation (modo reverso) es muy eficiente cuando tenemos muchas variables de entrada (parámetros) y una única salida escalar (la pérdida), ya que calcula todos los gradientes \(\frac{\partial L}{\partial \theta_i}\) con un coste computacional bajo (aproximadamente entre 2 y 3 veces el coste del forward pass). Por eso es el método estándar para entrenar redes neuronales.

# Función de activación y su derivada
def relu(x):
    return np.maximum(0, x)
def drelu(x):
    return (x > 0).astype(int)

# Dimensiones de la red
dim_entrada = 1
dim_oculta = 2
dim_salida = 1

# Ejemplo de entrada y objetivo
x = np.array([[0.5]])  # forma (1, 1)
y_real = np.array([[1.0]])  # forma (1, 1)

# Parámetros (inicialización aleatoria para demostración)
W1 = np.random.randn(dim_oculta, dim_entrada)  # (2, 1)
b1 = np.random.randn(dim_oculta, 1)            # (2, 1)
W2 = np.random.randn(dim_oculta, dim_oculta)   # (2, 2)
b2 = np.random.randn(dim_oculta, 1)            # (2, 1)
W3 = np.random.randn(dim_salida, dim_oculta)   # (1, 2)
b3 = np.random.randn(dim_salida, 1)            # (1, 1)

# ----- Paso hacia adelante (forward) -----
z1 = W1 @ x + b1         # (2, 1)
a1 = relu(z1)            # (2, 1)
z2 = W2 @ a1 + b2        # (2, 1)
a2 = relu(z2)            # (2, 1)
z3 = W3 @ a2 + b3        # (1, 1)
y_pred = z3              # (1, 1)

# Pérdida: Error cuadrático medio
perdida = 0.5 * np.sum((y_pred - y_real) ** 2)

# ----- Paso hacia atrás (backward) -----
# dJ/dy_pred
dJ_dy_pred = y_pred - y_real  # (1, 1)

# Gradientes de la capa de salida
# y_pred = z3
dJ_dz3 = dJ_dy_pred           # (1, 1)
dJ_dW3 = dJ_dz3 @ a2.T        # (1, 2)
dJ_db3 = dJ_dz3               # (1, 1)

# Segunda capa oculta
da2_dz2 = drelu(z2)           # (2, 1)
dJ_da2 = W3.T @ dJ_dz3        # (2, 1)
dJ_dz2 = dJ_da2 * da2_dz2     # (2, 1)
dJ_dW2 = dJ_dz2 @ a1.T        # (2, 2)
dJ_db2 = dJ_dz2               # (2, 1)

# Primera capa oculta
da1_dz1 = drelu(z1)           # (2, 1)
dJ_da1 = W2.T @ dJ_dz2        # (2, 1)
dJ_dz1 = dJ_da1 * da1_dz1     # (2, 1)
dJ_dW1 = dJ_dz1 @ x.T         # (2, 1)
dJ_db1 = dJ_dz1               # (2, 1)

# Ahora dJ_dW1, dJ_db1, dJ_dW2, dJ_db2, dJ_dW3, dJ_db3 son los gradientes para todos los parámetros
print("dJ/dW1:", dJ_dW1)
print("dJ/db1:", dJ_db1)
print("dJ/dW2:", dJ_dW2)
print("dJ/db2:", dJ_db2)
print("dJ/dW3:", dJ_dW3)
print("dJ/db3:", dJ_db3)

(Nota: Librerías modernas como PyTorch, TensorFlow y JAX implementan AD de forma automática. El usuario define el paso hacia adelante (la estructura del modelo y el cálculo), y la librería se encarga de calcular los gradientes eficientemente usando backpropagation.)

Inicialización de parámetros

La forma en que inicializamos los pesos (\(W\)) y sesgos (\(b\)) de una red neuronal es sorprendentemente importante, sobre todo en redes profundas (con muchas capas). Una mala inicialización puede impedir que la red aprenda.

No podemos inicializar todos los pesos y sesgos a cero, ya que permanecerán iguales a lo largo del entrenamiento. En cambio inicializamos todos los sesgos a cero y los pesos de forma aleatoria. Pero esto tiene una sutileza.

El problema de la varianza grande o pequeña

La clave está en controlar la varianza de las salidas de cada capa (activaciones) y de los gradientes que fluyen hacia atrás. Aquí, \(D_\ell\) representa el ‘fan-in’ de la capa \(\ell\), es decir, el número de neuronas en la capa anterior que alimentan a una neurona de la capa \(\ell\).

Los pesos se inicializan de forma aleatoria, lo que significa que los pesos tendrán una varianza \(\sigma^2\), que podemos escoger y puede ser distinta para cada capa. Resulta que el valor de esta varianza es importante para el entrenamiento de redes profundas ya que:

  • Varianza demasiado grande (e.g., \(\sigma^2 = 5 / D_\ell\)): Las activaciones tienden a crecer exponencialmente capa tras capa (activaciones explosivas). Como la activación de cada capa se multiplica por la de la capa anterior, esto puede llevar a la saturación de las funciones de activación (como sigmoid o tanh, o incluso valores enormes con ReLU) y a gradientes gigantescos (gradientes explosivos), desestabilizando el entrenamiento.
  • Varianza demasiado pequeña (e.g., \(\sigma^2 = 0.1 / D_\ell\)): Las activaciones tienden a atenuarse exponencialmente capa tras capa (activaciones desvanecientes). Como la activación de cada capa se multiplica por la de la capa anterior, esto provoca que las señales (y sus gradientes) que llegan a las primeras capas sean minúsculas (gradientes desvanecientes), impidiendo que esas capas aprendan (como el gradiente es pequeño, los pesos no se actualizan).

La @fig_varianza_activaciones ilustra cómo evoluciona la varianza de las activaciones a través de las capas para diferentes inicializaciones en una red con activación ReLU:

  • Pequeña (\(\sigma^2 = 0.1 / D_\ell\)): La varianza de las activaciones se atenúa rápidamente hacia cero.
  • Xavier/Glorot (\(\sigma^2 = 1 / D_\ell\)): Mantiene la varianza más estable, pero está diseñado teóricamente para activaciones lineales o simétricas (como tanh). Con ReLU, tiende a reducir la varianza. (La versión original de Xavier promedia fan-in y fan-out: \(\sigma^2 = 2 / (D_\ell + D_{\ell+1})\)).
  • He (\(\sigma^2 = 2 / D_\ell\)): Diseñada específicamente para ReLU. Mantiene la varianza de las activaciones aproximadamente constante a lo largo de las capas. El factor 2 extra (comparado con Xavier) compensa el hecho de que ReLU anula ~la mitad de las entradas (las negativas), lo que de otro modo reduciría la varianza.
  • Grande (\(\sigma^2 = 5 / D_\ell\)): La varianza explota brutalmente, saliéndose de escala y haciendo inviable el aprendizaje.
Código 
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# Parameters
n_layers = 20
fan_in = 100  # D_\ell
n_samples = 1000

# Variance schemes to compare
init_schemes = {
    r"Pequeña ($0.1/D_\ell$)": 0.1,
    r"Xavier/Glorot ($1/D_\ell$)": 1.0,
    r"He ($2/D_\ell$)": 2.0,
    r"Grande ($5/D_\ell$)": 5.0,
}
colors = ["gray", "blue", "green", "red"]

# Store results
all_variances = {}

# Simulate for each scheme
for (label, scale), color in zip(init_schemes.items(), colors):
    variances = []
    # Input: standard normal
    x = np.random.randn(fan_in, n_samples)
    var = np.var(x)
    variances.append(var)
    for layer in range(n_layers):
        # Initialize weights for this layer
        W = np.random.randn(fan_in, fan_in) * np.sqrt(scale / fan_in)
        # Linear + ReLU
        x = W @ x
        x = np.maximum(0, x)
        var = np.var(x)
        variances.append(var)
    all_variances[label] = variances

# Plot
plt.figure(figsize=(8, 5))
for label, color in zip(init_schemes.keys(), colors):
    plt.plot(all_variances[label], label=label, color=color, linewidth=2)
plt.xlabel("Capa (layer)")
plt.ylabel("Varianza de activaciones")
plt.title("Evolución de la varianza de activaciones según inicialización (ReLU)")
plt.legend()
plt.yscale('log')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Varianza de las activaciones según inicialización

Esto se puede entender porque la varianza de \(z^{(\ell)}\) está dada por:

\[ \sigma^2_{z^{(\ell)}} = \frac{1}{2}D_\ell \sigma_{\ell}^2\sigma_{z^{(\ell - 1)}}^2\,, \] {#eq_varianza_z}

donde \(\sigma_\ell^2\) es la varianza de los pesos de la capa \(\ell\). Esta relación también nos da una pista de cómo resolver el problema.

Queremos encontrar una relación entre la varianza de las preactivaciones de la capa \(\ell\) y la varianza de las activaciones de la capa \(\ell - 1\). Sea \(\sigma_{\ell}^2\) la varianza de los pesos de la capa \(\ell\) y \(\sigma_{z^{(\ell)}}^2\) la varianza de las preactivaciones de la capa \(\ell\).

Para empezar, calculamos el valor esperado de la preactivación \(z^{(\ell)}\) en la capa \(\ell\):

\[ \langle z^{(\ell)} \rangle = \langle W^{(\ell)} \cdot a^{(\ell-1)} + b^{(\ell)} \rangle = \langle W^{(\ell)} \rangle \cdot \langle a^{(\ell-1)} \rangle = 0\,. \]

Entonces su varianza es:

\[ \text{Var}(z_{i}^{(\ell)}) = \langle (z_{i}^{(\ell)})^2 \rangle = \left\langle \left(\sum_{j=1}^{D_\ell} W_{ij}^{(\ell)} a_{j}^{(\ell-1)}\right)^2 \right\rangle = \sum_{j = 1}^{D_\ell} \langle (W_{ij}^{(\ell)})^2 \rangle \langle (a_{j}^{(\ell-1)})^2 \rangle = \sigma_{\ell}^2 \sum_{j=1}^{D_\ell}\langle (a_{j}^{(\ell-1)})^2 \rangle\,. \]

Ahora bien, como \(ReLU\) corta a cero las entradas negativas, la varianza de las activaciones \(a^{(\ell)}\) es la mitad de la varianza de \(z^{(\ell)}\) (ya que la integral de \(x^2\) desde 0 hasta \(\infty\) es la mitad de la integral desde \(-\infty\) hasta \(\infty\)). Entonces:

\[ \sigma_{z^{(\ell)}}^2 = \frac{1}{2}D_\ell \sigma_{\ell}^2\sigma_{z^{(\ell - 1)}}^2\,. \]

Inicialización para el forward pass: \(\sigma_\ell^2 = 2/D_\ell\)

Para asegurar que la varianza de las activaciones \(a^{(\ell)}\) se mantenga aproximadamente constante durante el forward pass (al pasar de la capa \(\ell\) a \(\ell+1\)) cuando se usa la activación ReLU, Kaiming He et al. propusieron inicializar los pesos \(W^{(\ell)}\) muestreando de una distribución con media 0 y varianza:

\[ \sigma_\ell^2 = \text{Var}(W_{ij}^{(\ell)}) = \frac{2}{D_\ell} \]

donde \(D_\ell\) es el fan-in de la capa \(\ell\). Si reemplazamos esto en la @eq_varianza_z obtenemos \(\sigma_{z^{(\ell)}}^2 = \sigma_{z^{(\ell - 1)}}^2\), tal que la varianza es estable.

Inicialización para el backward pass: \(\sigma^2 = 2/D_{\ell+1}\)

De forma análoga, si analizamos el backward pass, para que la varianza de los gradientes \(\frac{\partial \mathcal{J}}{\partial z^{(\ell)}}\) se mantenga estable al retropropagar desde la capa \(\ell+1\) a la capa \(\ell\) (asumiendo activación ReLU), la condición requiere una varianza basada en el ‘fan-out’:

\[ \sigma^2 = \text{Var}(W_{ij}^{(\ell)}) = \frac{2}{D_{\ell+1}} \]

Notamos que la condición ideal para el forward pass (usa \(D_h\)) y para el backward pass (usa \(D_{h+1}\)) solo coinciden si \(D_h = D_{h+1}\). En la práctica, se puede usar la inicialización de He cuando la activación es ReLU, o un promedio \(\sigma_\ell^2 = \frac{4}{D_h + D_{h+1}}\).

Nota: Cuando la activación es tanh o sigmoide, la derivación de la relación entre las varianzas sugiere usar \(\sigma_\ell^2 = \frac{1}{D_h}\), que se llama inicialización de Xavier/Glorot.

Ejemplo de código de entrenamiento

Existen librerías que implementan todos los pasos de este entrenamiento. La más popular es tal vez PyTorch. Este curso no está enfocado en implementar el código ya que este cambia constantemente, y el simple uso de librerías se puede lograr leyendo la documentación, pidiendo ayuda en foros como StackOverflow o preguntándole a una IA. Sin embargo, damos algunos ejemplos, empezando por un ejemplo de cómo se implementaría lo visto hasta ahora usando PyTorch.

import torch, torch.nn as nn
from torch.utils.data import TensorDataset, DataLoader
from torch.optim.lr_scheduler import StepLR

# Definir el tamaño de la entrada, el tamaño de las capas ocultas y el tamaño de la salida
D_i, D_k, D_o = 10, 40, 5
# crear el modelo con dos capas ocultas
model = nn.Sequential(
    nn.Linear(D_i, D_k),
    nn.ReLU(),
    nn.Linear(D_k, D_k),
    nn.ReLU(),
    nn.Linear(D_k, D_o)
)

# Inicialización de He
def ini_pesos(capa_in):
    if isinstance(capa_in, nn.Linear):
        nn.init.kaiming_normal_(capa_in.weight, nonlinearity='relu')
        capa_in.bias.data.fill_(0.0)
model.apply(ini_pesos)

# Función de pérdida ECM
criterion = nn.MSELoss()
# Construir el optimizador por SDG e inicializar la tasa de aprendizaje y el momentum
optimizer = torch.optim.SGD(model.parameters(), lr=0.1, momentum=0.9)
# Objeto que reduce la tasa de aprendizaje por la mitad cada 10 épocas
scheduler = StepLR(optimizer, step_size=10, gamma=0.5)

# Crear 100 muestras aleatorias
x = torch.randn(100, D_i)
y = torch.randn(100, D_o)
dataloader = DataLoader(TensorDataset(x, y), batch_size=10, shuffle=True)

# Entrenar el modelo, pasando 100 veces por los datos
for epoch in range(100):
    epoch_loss = 0.0
    for i, data in enumerate(dataloader):
        x_batch, y_batch = data
        # Poner a cero los gradientes
        optimizer.zero_grad()
        # Forward pass
        y_pred = model(x_batch)
        # Calcular la pérdida
        loss = criterion(y_pred, y_batch)
        # Backward pass
        loss.backward()
        # Actualizar los pesos
        optimizer.step()
        # Actualizar la pérdida
        epoch_loss += loss.item()
    # Imprimir la pérdida
    print(f"Epoch {epoch:5d}, Loss: {epoch_loss:.3f}")
    # Actualizar la tasa de aprendizaje
    scheduler.step()

Ejemplo aplicado de una red neuronal

En el libro hay un ejemplo para un problema de clasificación con 40 características. Recomiendo leerlo.

Otro ejemplo típico es la clasificación de dígitos escritos, el llamado MNIST. Consiste en un conjunto de imágenes en blanco y negro que representan dígitos escritos a mano. Se entrena una red para identificar a cuál dígito corresponde cada imagen.

Primero cargamos los datos. Note que necesitamos normalizarlos. En general las redes neuronales funcionan mejor si los datos tienen media cero y varianza unitaria. Por eso dividimos la media y la desviación estándar del conjunto de entrenamiento y usamos esas para normalizar el conjunto de entrenamiento y el de prueba.

Para lograrlo, primero cargamos los datos sin normalizar, les calculamos su media y varianza, y finalmente usamos transforms.Normalize para normalizar los datos.

import torch
import torch.nn as nn
import torch.optim as optim
from torchvision import datasets, transforms
import matplotlib.pyplot as plt

# Definir transformaciones para los datos
transform = transforms.Compose([
    transforms.ToTensor()
])

# Cargar datos de MNIST
train_dataset = datasets.MNIST('./data', train=True, download=True, transform=transform)
test_dataset = datasets.MNIST('./data', train=False, transform=transform)

# Calcular media y desviación estándar del conjunto de entrenamiento
train_loader_temp = torch.utils.data.DataLoader(train_dataset, batch_size=len(train_dataset))
data = next(iter(train_loader_temp))[0]
mean = data.mean().item()
std = data.std().item()

# Aplicar normalización con los valores calculados
transform_norm = transforms.Compose([
    transforms.ToTensor(),
    transforms.Normalize((mean,), (std,))
])

# Recargar los datasets con la normalización correcta
train_dataset = datasets.MNIST('./data', train=True, download=True, transform=transform_norm)
test_dataset = datasets.MNIST('./data', train=False, transform=transform_norm)

# Definir dataloaders
train_loader = torch.utils.data.DataLoader(train_dataset, batch_size=64, shuffle=True)
test_loader = torch.utils.data.DataLoader(test_dataset, batch_size=1000, shuffle=False)

Podemos graficar algunas de las imágenes del conjunto de entrenamiento.

# Graficar algunas imágenes del conjunto de entrenamiento
import matplotlib.pyplot as plt

# Tomar 16 imágenes del conjunto de entrenamiento
images, labels = next(iter(train_loader))
images = images[:16]
labels = labels[:16]

# Graficar las imágenes
plt.figure(figsize=(10, 10))
for i in range(16):
    plt.subplot(4, 4, i + 1)
    plt.imshow(images[i].numpy().squeeze(), cmap='gray')
    plt.title(labels[i].item())
    plt.axis('off')
plt.show()

Ahora vamos a entrenar un perceptrón multicapa para clasificar las imágenes.

# Definir una red neuronal MLP simple
class MLP(nn.Module):
    def __init__(self):
        super(MLP, self).__init__()
        self.flatten = nn.Flatten()
        self.layers = nn.Sequential(
            nn.Linear(28 * 28, 128),
            nn.ReLU(),
            nn.Linear(128, 64),
            nn.ReLU(),
            nn.Linear(64, 10)
        )
    
    def forward(self, x):
        x = self.flatten(x)
        return self.layers(x)

# Inicializar modelo, función de pérdida y optimizador
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
# # Diagnóstico para disponibilidad de la GPU
# print("Calculando en: ", device)
# print("Dispositivo disponible: ", torch.cuda.is_available())
# if torch.cuda.is_available():
#     print("Dispositivo actual: ", torch.cuda.current_device())
#     print("Nombre del dispositivo: ", torch.cuda.get_device_name(torch.cuda.current_device()))
model = MLP().to(device)
criterion = nn.CrossEntropyLoss()
optimizer = optim.SGD(model.parameters(), lr=0.01, momentum=0.9)

# Función para evaluar el modelo
def evaluate(model, data_loader, criterion, device):
    model.eval()
    total_loss = 0
    correct = 0
    total = 0
    
    with torch.no_grad():
        for data, target in data_loader:
            data, target = data.to(device), target.to(device)
            output = model(data)
            total_loss += criterion(output, target).item() * data.size(0)
            
            _, predicted = torch.max(output.data, 1)
            total += target.size(0)
            correct += (predicted == target).sum().item()
    
    return total_loss / total, correct / total

# Entrenamiento
epochs = 10
train_losses = []
test_losses = []
train_accuracies = []
test_accuracies = []

for epoch in range(epochs):
    model.train()
    running_loss = 0.0
    correct = 0
    total = 0
    
    for data, target in train_loader:
        data, target = data.to(device), target.to(device)

        optimizer.zero_grad()
        output = model(data)
        loss = criterion(output, target)
        loss.backward()
        optimizer.step()
        
        running_loss += loss.item() * data.size(0)
        _, predicted = torch.max(output.data, 1)
        total += target.size(0)
        correct += (predicted == target).sum().item()
    
    train_loss = running_loss / len(train_loader.dataset)
    train_accuracy = correct / total
    
    # Evaluar en conjunto de prueba
    test_loss, test_accuracy = evaluate(model, test_loader, criterion, device)
    
    # Guardar métricas
    train_losses.append(train_loss)
    test_losses.append(test_loss)
    train_accuracies.append(train_accuracy)
    test_accuracies.append(test_accuracy)
    
    # # Para ver el progreso del entrenamiento
    # # comentarlo a la hora de producir el sitio web
    # print(f'Epoch {epoch+1}/{epochs}: '
    #       f'Train Loss: {train_loss:.4f}, Train Acc: {train_accuracy:.4f}, '
    #       f'Test Loss: {test_loss:.4f}, Test Acc: {test_accuracy:.4f}')

Después de entrenar, graficamos las pérdidas y la precisión en el conjunto de entrenamiento y prueba.

# Graficar resultados
plt.figure(figsize=(10, 6))

plt.subplot(1, 2, 1)
plt.plot(train_losses, label='Pérdida en entrenamiento')
plt.plot(test_losses, label='Pérdida en prueba')
plt.xlabel('Época')
plt.ylabel('Pérdida')
plt.legend()
plt.title('Pérdida vs. Época')

plt.subplot(1, 2, 2)
plt.plot(train_accuracies, label='Precisión en entrenamiento')
plt.plot(test_accuracies, label='Precisión en prueba')
plt.xlabel('Época')
plt.ylabel('Precisión')
plt.legend()
plt.title('Precisión vs. Época')

plt.tight_layout()
plt.show()

Ejercicios Sugeridos

6.4, 6.8, 6.9, 7.7, 7.11, 7.12